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FOREWORD 


The  work  described  in  this  report  was  performed  in  the  Fire  Control 
Formulation  Branch  (K41),  SLBM  Research  and  Analysis  Division,  Strategic- 
Systems  Department,  and  was  authorized  under  Strategic  Systems  Program  Office 
Task  Assignment  36401.  This  work  was  necessitated  by  the  need  to  formulate  an 
exact  method  of  computing  geodetic  latitude  and  altitude. 

This  report  has  been  reviewed  and  approved  by  Davis  Owen;  Johnny  Boyles; 
Robert  Gates,  Head,  Fire  Control  Formulation  Branch;  and  Sheila  Young,  Head, 
SLBM  Research  and  Analysis  Division. 

This  document  is  a  revised  version  of  NSWC  TR  85-85,  dated  April  1985.  The 
distribution  statement  has  been  changed. 

Approved  by; 


r' —  .  •  '77 ' 

R.  L.  SCHMIDT,  Head 

Strategic  and  Space  Systems  Department 
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INTRODUCTION 


Fallowing  is  a  simple  and  efficient  model  for  calculating  the  exact 
geodetic  latitude  and  altitude  of  an  arbitrary  point  in  space,  given  the 
coordinates  of  that  point.  The  model  assumes  the  earth  to  be  an  oblate 
spheroid;  i.e.,  an  ellipsoid  of  revolution  whose  semimajor  axis  (a)  is  the 
radius  of  the  circle  described  by  the  equatorial  plane,  and  whose  semituinot 
axis  (b)  is  a  line  joining  its  center  and  one  of  its  poles.  The  point  in 
question  can  be  either  inside  or  outside  the  ellipsoid,  but  not  too  close  to 
the  center  of  the  ellipsoid.  The  solution  of  the  problem  and  its  constraints 
will  now  be  expounded. 


FORMULATION 


Choose  an  earth-fixed  Cartesian  coordinate  system  whose  origin  coincides 
with  the  center  of  the  ellipsoid.  The  unit  vectors  i,  j,  k  coincide 
with  the  (x,  v,  z)  axes,  respectively.  The  +  z  axis  points  in  the  direction  of 
the  North  Pole.  The  +x  axis  is  the  line  of  intersection  of  the  equatorial 
plane  and  the  plane  of  zero  longitude.  The  +y  axis  completes  a  right-handed 
coordinate  system. 

An  equation  of  the  ellipsoid  in  this  frame  is 

~-r  +  m  1 ,  where  R  *  V  x2  +  y2  .  (  1  ) 

a4  b4 

Let  P(x0,  yQ,  z0)  be  the  coordinates  of  the  given  point.  It  is  desired  to 
find  the  points  P(x,  y,  z)  at  which  the  normals  from  P  to  the  surface  of  the 
ellipsoid  cut  the  ellipsoid. 


The  slope  of  a  normal  to  the  ellipsoid  at  any  point  on  its  surface  is 
given  by 


1  ^  a2  z 

dz  b2  R 
dR 


d2 

where  —  +  r-r  =  1  . 
a4  b4 


(2) 


Therefore  the  slope  of  a  normal  from  P  is 


z0  -  z 
Ro  -  R 


•7 


a~  z 
b2  R 


where  R^, 


v: 


+  V4 
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i-»-  , 


(R^  -  R)a2z  *  (  z0  -  z)b2R  , 


Rfl  a2  z  =  R{a2z  +  b2  (  z0  -  z)}  =  R{(a2  -  b2)z  +  b2  z0  J  . 

Squaring  both  sides  and  expressing  R  in  terms  of  z, 
a2b2R^z2  *  (b2  -  z2){(a“  -  b2)z  +  b2  z0  J  . 

Writing  the  above  equation  in  descending  powers  of  z,  we  obtain  the 
following  biquadratic  in  z: 

(a2  -  b2  )  z*  +  2b2(a2  -  b2)z0z3  +  b2{a2R^  +  b2zg  -  (a2  -  b2)2)z2  - 
2 b**  (a2  -  b2)z02  -  b6  z|  *  0. 

Now  z  is  either  positive  or  negative  and  therefore  is  not  an  appropriate 
variable  for  the  solution.  However,  this  difficulty  is  easily  circumvented  by 

z 

introducing  a  new  dimensionless  variable  k  ■  — .  Using  the  constraint  that  z 

zo 

always  has  the  same  sign  as  z0  ,  we  see  that  k  is  always  positive.  Putting  z  » 
z0 k  in  the  biquadratic  and  writing  the  resulting  equation  such  that  the 
coefficient  of  k1*  is  unity,  we  finally  obtain 


k1*  +  2pk3  +  qk2  +  2rk  +  s  *  0 


where 


b2 

p  -  ^r-r-bT 

b2{a2R^  +  b2z^  -  fa2  -  b2)2} 
q  c  (a2  -  b2)2z2 


(a2  -  b2 )zl 

b6 _ 

(a2  -  b2  ) 2  Zq 


Now  we  will  emplov  two  of  thrpo  powerful  theorems  in  the  Theorv  of  F.qua- 
tions  to  expose  the  nature  of  the  roots  of  (*).  The  third  theorem  will  be 
needed  later  on.  The  proofs  of  these  theorems  are  given  in  Reference  1. 
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I.  An  equation  f(x)  =  0  cannot  have  more  positive  roots  than  there  are 

changes  of  sign  in  f(x),  and  cannot  have  more  negative  roots  than  there  are 

changes  of  sign  in  f(-x). 

II.  Every  equation  which  is  of  an  even  degree  and  has  its  last  term 
negative  has  at  least  two  real  roots,  one  positive  and  one  negative. 

III.  Every  equation  of  an  odd  degree  has  at  least  one  real  root  whose  sign 

is  opposite  to  that  of  its  last  term- 

By  observing  (*)  and  (3),  (4),  (5),  (6),  we  deduce  at  onre  the  following 
conclusions : 

(i)  Since  the  last  term  of  (*),  namely  s,  is  negative,  there  are  at  least 
two  real  roots  of  opposite  sign  (using  IT). 

(ii)  Since  there  is  only  one  change  of  sign  in  (*),  irrespective  of 
whether  q  is  positive  or  negative,  there  is  at  most  one  positive  root 
(using  I). 

From  (i)  and  (ii)  we  see  immediately  that  (*'  has  exactly  one  positive 
root.  Since  k  is  constrained  to  be  positive,  this  positive  root  is  the  one  we 
seek . 


SOLUTION  OF  BIQUADRATIC 

The  solution  of  (*)  is  effected  by  a  standard  method  known  as  Ferrari 's 
method  which  will  now  be  enunciated. 

In  the  equation 

k1*  +  2pk3  +  qk2  +  2rk  +  s  *  0 

add  to  each  side  (ck  +  d)2,  the  quantities  c  and  d  being  determined  so  as  to 
make  the  left-hand  side  a  perfect  square;  then 

k1*  +  2pk3  +  (q  +  c2  )k2  +  2(r  +  cd)k  +  s  +  d2  *  ( ck  +  d)2  . 

Suppose  that  the  left  side  of  the  equation  is  equal  to  ( k2  +  pk  +  t)2  , 
then  by  comparing  the  coefficients,  we  have 

p2  +  2t  =  q  +  c2 ,  pt  *  r  +  cd,  t2  =  s  +  d2 .  (71 

By  eliminating  c  and  d  from  these  equations,  we  obtain 

( pt  -  r)2  =  (2t  +  pJ  -  q)(t2  -  s), 

or  2t3  -  qt2  +  2( pr  -  s)t  -  p2 s  +  qs  -  r2  =0.  (8) 
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From  (3),  (4),  (5),  and  (6),  we  see  that 
pr  -  s  =  0,  -p2 s  +  qs  -  r2  = 


a2  b8  R^j 


(a2  -  b2  T  zl 


Substituting  into  (8),  we  obtain  the  following  cubic  in  t: 
a2b5 

2t3  -  qt2  -  - — Tt  ~ — —  =  0  . 


(a2  -  b2  )' 


(  8’  ) 


Applying  Theorem  III  to  (8'),  we  see  that  (8')  has  at  least  one  positive 
root.  Now  applying  Theorem  I  to  (8'),  we  see  that  it  has  at  most  one  positive 
root,  irrespective  of  whether  q  is  positive  or  negative.  Hence  :8')  has 
exactly  one  positive  root. 

The  solution  of  a  cubic  is  accomplished  by  a  standard  method  known  as 
Cardan's  Solution. 

First  eliminate  the  t2  term  in  (8')  by  making  the  substitution 


t  -  t' 


i.e. 


> 


2(t'  +  |)3  _  q(t'  +  a)2 


a2b8R2 

(a2  -  b2 11*  zl  ® 


or 


t'3 


_  aL.  . 

108 


a2b8Ro2 

2(  a2  -  b2  )h  zl 


0  . 


Let 


f  = 


,  h 


.  aL.  _ 

108 


a2  b8  R2 

27  a 2  -  b2  r  zl 


so  that  the  above  equation  can  be  written  in  the  form 

t ' 3  +  f  t '  +  h  *  0  .  (9) 

To  solve  (9),  lot  t'  -  u  +  v;  then 

t>3  =  < u  +  v ) 3  *  u3  +  v3  +  3uv  (u  +  v)  *  u3  +  v3  +  3uvt' 
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and  ( 9 )  becomes 

uJ  +  v3  +  (  3uv  +f)t'+h=0. 


At  present,  u  and  v  are  any  two  quantities  subject  to  the  condition  that 
their  sum  is  equal  to  one  of  the  roots  of  (9);  if  we  further  suppose  that  thev 
satisfy  the  equation  3uv  +  f  *  0,  they  are  completely  determinate.  We  thus 
obtain 


f  3 

u3  +  v3  =  -  h,  u3v3  *  -  — 

Hence,  u3  ,  v3  are  the  roots  of  the  quadratic 


a 2  +  ha  -  ttz  =  0 


Thus , 


h 

°  =  ~  2 

f3 

'  +  27 

Putting  u3  «= 

-?♦ 

/  Hi  +  fi 

V  4  27 

t 

1 

fi 

> 

/f7! 

the  relation 

t '  =  u 

+  V. 

Hence,  t '  * 

(-!* 

/Hi  *  it 

y  4  27 

1*  3  /  h 

)  +  0  2  -  , 

/F7!) 

h2  f3 

The  solution  (10)  valid  provided  that  - —  +  i  0.  It  will  now  be 

shown  that  this  constraint  is  valid  for  all  points  of  interest. 
h2  fi  l  i  ni  a2b8R^  ,2  6 


hi  fi  I  (  q3  *  u  ^ 

4  27  *4)  108  2(a2  -  b2  )" 

a2  b8  R|  q3  i 


i2 _ a 

)  46656 


a2  b8  Rj  q3  a’  b1  6  Rj 

432(a2  -  b2  )h  Zp  +  16(a2  -  b2  )8  zjj 

a2  b1  H  { a2  Rg  +  b2  z|  -  (a2  -  b2)2}3 


432(a2  -  b2  )!0z30 


a1*  b1  6  Rg 

T6(a2  “  bM8'zfr  uslng  <M 


a2  b1 P.q 


432( a2 


_ ‘  ~j  i  o  7_ i  o  i  Ro  +  h2zo  ’  (a?  "  b:>:P  + 


27a2b2(a2  -  b2)2Rg’z/jl 


NSWf  TH 


Thus,  the  required  condition  is 

{a-R2  +  b2z2  -  (a2  -  b2  +  2  7 a 2  b2  (  a 2  -  b2l2^ij  >  0  .  Mi' 

Bv  inspection  and  also  by  trial,  it  is  seen  that  the  above  <-onstraint  is 
valid  at  all  points  of  interest.  As  a  matter  of  fart,  the  only  points  where  it 
does  not  hold  are  those  that  are  relatively  close  to  the  origin,  i  ,  the 
center  of  the  ellipsoid;  and  these  points  are  of  no  practical  interest. 


a2  b5  R5 

2 (a2  -  h2  )7zl 


bHa-’R2  +  b2z2  -  'a2  -  b-  )2p 
1 08 (  a2  -  b2  )b  'zjj 


^R2 

2< a2  -  b2  >’  zl 


108( a2  -  b2  )fa  z □ 


[  {  a2  Rq  +  b2  Zq 


a2  -  b2  )2  P  +  54a2  b2 1  a2  -  hi:)2R^z^ 


If  the  constraint  (11)  holds,  then  the  expression  in  the  brackets  above  is 
positive.  Hence,  h  is  negative.  Applying  Theorem  III  to  (  9 ) ,  we  conclude  that 
(9)  has  at  least  one  positive  root.  Now  applying  Theorem  I  to  (9)  and  using 
the  fact  that  both  f  and  h  are  negative,  we  conclude  that  (9)  has  at  most  one 
positive  root.  Hence,  (9)  has  exactLv  one  positive  root.  It  will  now  be  shown 
that  this  positive  root  is  given  by  (10). 

_  h2  f3 

Let  g  *  +  T7  ♦  then 

u.(vf-|)‘/3.  v.(-  Vf-|)‘/3  • 

u2  f 3  q2 

Now  g  -  —  *  —  <  0,  since  f  =»  -  -pj  <  0 

Therefore.  (  VF  -  <  V~g  +  -)  <  0,  which  implies  that  VT  -  -  and  Vg"  +  “ 

are  of  opposite  sign;  however,  since  h  <  0,  Vg"  -  cj  >  Vg  +  7  -  Hence,  Vg"  -  ^  > 


0  >  Vg  +  — ;  i.e.,  u  >  v  >  0.  Thus,  t'  -  u  +  v  >  0. 


y.Wf 


Ik 

Proceeding  with  our  derivation,  w*>  have 

t  -  t'  +  a  >  0  . 

b 

Solving  foi  c  and  d  from  the  system  of  equations  (  7),  we  have 


sow  (  k2  +  pk  +  t)-  =  tck  +  d>-  , 


or  k2  +  pk  +  t  =  *  ( ok  +  d ! 
from  which  we  obtain  the  two  quadratics  in  k: 

k'+(p-c)k  +  t-da0  (  !  2  ' 
k2+(p  +  c)k  +  t+  d  =  0  (lii 
From  (12)  and  (13),  we  ohtain  the  following  four  roots  of  <*): 
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It  has  already  been  shown  that  (*)  has  exactly  one  positive  root.  Since 
t.2  -  d2  =  s  <  0,  <t  +  d)(t.  -  d)  <  0,  which  implies  that  t  +  d  and  t.  -  d  are  of 
opposite  sign;  blit  d  >  0;  thus  t+d>()>t-d.  Now  t  -  d  i?  the  last  term 
of  (12)  and  has  just  been  shown  to  be  negative.  Hence,  bv  applying  Theorem  II 
to  (12),  we  see  immediately  that  (12)  must  contain  the  desired  positive  root. 
Inspecting  the  first  two  roots  above,  it  is  clear  that  k,  is  the  required 
positive  root. 

The  rest  nf  the  solution  follows  trivially.  Let  k  =  k,  ,  then 
7  =  70k  , 


KSWC  TK  8^-8 S 


R  *  a  / 1  -  rr  ,  and 
b" 


x  -  xA 


R« 


y  “  y°  r”  ’  r  wherP  Ro  *  o 


\  *  tan  1  *- 
x 


J 


(Here  \  is  the  longitude  at  (x,  yl). 

The  above  solution  is  not  applicable  to  the  case  ?0  =  0  since  two  of  the 
bacic  parameters  of  the  solution,  q  and  s,  involve  a  division  by  z0  .  Also,  the 
sol  Jtion  is  not  applicable  to  the  case  R0  «  0.  Hence  these  two  cases  have  to 
be  created  separately;  nevertheless,  this  poses  no  difficulty  since  the  results 
are  already  known  then. 


The  geodetic  altitude  is  given  by 


.  (  R0  Z0  \  /  2 

H  =  sign  I  ~  ~  ~  1  )  V  (x0  -  x)z  + 


%~>'r  +  (zo 


z)‘ 


=  sign  (  ~°  +  ^  -  1 ) 

a  b  ' 

In  order  to  find  the  geodetic  latitude  and  the  unit  normal  vector  at  (x, 
y,  t),  it  is  desirable  to  introduce  a  geometric  term,  Ne,  which  is  never  zero. 
Ne  is  defined  to  be  the  distance  along  the  ellipsoidal  normal  from  the  surface 
of  the  ellipsoid  to  the  z-axis  (Reference  2).  (See  Figure  1.) 


FIGURE  1.  ELLIPSOIDAL  NORMAL 
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From  Figure  1 ,  we  have 
R 


COS  0 


Ne 


a2  z 


Also , 

tan  0  a  prji  »  from  (2). 
Hence , 


Ne  =  R  sec  0  *  R 


^  1  +  tan‘0  *  R  >  1  +  ~ 


b"  R2 


-2  V  a'  z2  ±  h-  R2 


b2 


sin  0  =  cos  0  tan  o  *  rrr 

b-Ne 

Thus,  the  components  of  the  unit  normal  vector  at  (x,  y,  z)  are 

2 


H 


1  Ne 


’ H*  •  k  1 H3  ’  (b)  fe  • 


The  geodetic  latitude  0,  is 
0  =  sin"1  H3  . 

This  completes  our  solution. 


CONCLUSION 


In  contrast  to  other  methods  currently  in  use,  this  method  theoretically 
computes  the  exact  values  of  the  geodetic  latitude  and  altitude,  assuming  the 
earth  to  be  a  perfect  ellipsoid  of  revolution.  The  method  was  programmed  and 
run  for  many  given  values  of  x0 ,  yg ,  z0 .  It  gave  results  which  are  consistent 
with  those  of  other  methods,  including  Brookshire's  approximation  (Reference  3) 
and  an  iterative  technique  (Reference  4).  Roundoff  errors  are  negligible.  For 
all  practical  purposes,  the  results  are  excellent  and  the  method  is  definitely 
recommended  as  an  alternative  approach  to  problems  involving  geodetic  latitude. 
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APPENDIX 


THE  INVERSE  PROBLEM 


An  arbitrary  point  in  space  is  defined  by  A,  ©.  and  H.  Solve  for  (H, 
H2  ,  H3  ) ,  (x,  v,  z),  and  ( x0  ,  v0  ,  z0  ) . 

This  problem  is  relatively  easy  to  solve.  From  page  9,  we  have 

b2  . 

R  «  he  cos  ©,  z  *  — r  he  sin  ©. 

a~ 

R2  z2 

Using  the  relation  pr  +  p-  =  1 >  we  obtain 


Ne2  / 

a2  V 


cos2©  +  -y  sin2 


©^  ■  1  ; 


/  77  b2  .  . 

'cos2©  +  — y  am2' 
ac 


1  -  ( 1  -  ~r)  sin2 


Also,  from  page  8,  we  have 


X  »  tan-1  ^  ■  sin'1  ^ 


Therefore , 


x  *  R  cos  X  *  Ne  cos  ©  cos  X, 


y  ■  R  sin  A  *  Ne  cos  ©  sin  A. 


Hence , 


Ne  cos  ©  cos  A 
Ne  cos  ©  sin  A  , 

b2  K  .  . 

—  Ne  sin  © 


A-l 


.vsur  tr 


- 

X 

Ne 

COS  0  cos 

\’ 

h2 

sa 

y 

— 

cos  0  sin 

X 

Ne 

H,. 

.‘5 

■  2_ 

sin  o 

L. 

- 

Ne_ 

From  Figure  1,  it  follows  immediately  that 


xo 

X 

V 

( Ne  +  H)  cos  O  cos  X 

y0 

=S 

V 

+  H 

H, 

= 

( Ne  +  H)  cos  0  sin  X 
,  ■) 

_z0. 

z 

.Hj. 

(~  Ne  +  H)  sin  0 

A- 2 
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